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Reflection and refraction of electromagnetic waves by artificial periodic composites (metamateri- 
als) can be accurately modeled by an effective medium theory only if the boundary of the medium is 
explicitly taken into account and the two effective parameters of the medium - the index of refraction 
and the impedance - are correctly determined. Theories that consider infinite periodic composites 
do not satisfy the above condition. As a result, they cannot model reflection and transmission by 
finite samples with the desired accuracy and are not useful for design of metamaterial-based devices. 
As an instructive case in point, we consider the "current-driven" homogenization theory, which has 
recently gained popularity. We apply this theory to the case of one-dimensional periodic medium 
wherein both exact and homogenization results can be obtained analytically in closed form. We 
show that, beyond the well-understood zero-cell limit, the current-driven homogenization result is 
inconsistent with the exact refiection and transmission characteristics of the slab. 



I. INTRODUCTION 

In the past decade, interest in electromagnetic homog- 
enization theories has experienced a remarkable revival, 
especially when applied to artificial periodic composites 
(metamaterials)^^. The ultimate goal of any homoge- 
nization or effective medium theory (EMT) is to describe 
reflection and refraction of waves by finite samples. In the 
case of homogeneous natural materials, an accurate de- 
scription of this kind is possible only if both the index of 
refraction and the impedance of the material are known 
with sufficient precision. Correspondingly, all EMTs at- 
tempt to replace a composite sample with a sample of 
the same overall shape but spatially-uniform effective re- 
fractive index and impedance. However, when the EMTs 
are tested or evaluated, the attention is frequently paid 
only to the physical quantities that depend on the in- 
dex of refraction alone but not on the impedance. In 
particular, this is the case for all EMTs that consider 
infinite composites and do not account for the bound- 
ary of the medium. Still, these theories always predict 
some impedance, and the question remains whether this 
prediction is applicable to finite samples. 

The analysis is relatively simple in the classical ho- 
mogenization limit /i — >■ 0, where h is the heterogeneity 
scale such as the lattice period of a composite. Note that 
here we assume that all physical characteristics of the 
constituents of the composite arc independent of h. Wc 
will refer to this kind of EMT as "standard" . Note that 
an alternative approach has been proposed^i^ in which 
the limit /i — > is also taken but the permittivity of 
one of the composite constituents is assumed to depend 
on h. This theory is of a more general or, as we shall 
say, of the "extended" type. The fundamental differ- 
ences between standard and extended theories have been 
discussed by Bohren^. What is important here is that 



standard EMTs do not mix the electric and magnetic 
properties of the composite constituents^. This means, 
in particular, that the effective permeability obtained in 
a standard EMT is identically equal to unity if the con- 
stituents of the composite are intrinsically nonmagnetic. 
A closely related point is that, in standard theories, the 
impedance of the medium can be inferred from the bulk 
behavior of waves as long as we accept that the effective 
permeability is trivial. It can be proved independently 
that, in the /i — >■ limit, this choice of impedance is con- 
sistent with the exact Frcsncl reffcction and refraction 
coefficients at a planar boundary^. Thus, in a standard 
theory, both the impedance and the refractive index are 
consistent with reflection and refraction properties of a 
finite sample. 

However, standard EMTs are typically viewed as inad- 
equate in the modern research of electromagnetic meta- 
materials because these theories do not predict or de- 
scribe the phenomenon of "artificial magnetism" , which 
has a number of potentially groundbreaking applica- 
tions^. This difficulty is not characteristic of the ex- 
tended theories. An extended EMT either docs not cm- 
ploy the limit /i — > or, otherwise, assumes mathemat- 
ical dependence between h and other physical parame- 
ters of the composite. The main question we consider 
in this paper is whether an extended EMT can predict 
the refractive index and impedance simultaneously and 
in a reasonable way. Of course, a refractive index per 
se (generally, tensorial and dependent on the direction 
of the Bloch wave vector) can always be formally intro- 
duced for a Bloch wave. This can be done even in the 
case when the composite is obviously not electromagnet- 
ically homogeneous. But all extended EMTs yield both a 
refractive index and an impedance, and in the case of infi- 
nite unbounded media there is no way to tell whether this 
homogenization result is reasonable. In this paper, we 



present a case study by comparing the so-called current- 
driven homogenization theory (which is of extended type 
and is formulated for an infinite medium) to exact results 
in a layered finite slab. Note that, although we analyze 
a particular EMT, the central theme of this paper is re- 
lated to the fundamental difference between standard and 
extended EMTs. 

There are, of course, many extended EMTs currently 
in circulation. Theories of this kind have been first pro- 
posed by LewiniS and Khizhnyakii"— but they came 
to the fore more recently in the work of Niklasson et 
ali^, Doylei^, and Waterman and Pederseni^, who have 
generalized the classical Maxwell-Garnett approximation 
to account for the magnetic dipole moments of spher- 
ical particles (e.g., computed using Mie theory). Al- 
though the extended Maxwell-Garnett approximation of 
Refs. (Tj-lla applies only to the dilute case, it has served 
as an important precursor of several more generally ap- 
plicable extended EMTs. Among these we can men- 
tion the modified multiscale approach^, Bloch analysis 
of electromagnetic latticesii"— , coarse-graining (averag- 
ing) of the electromagnetic fields using curl-conforming 
and div-conforming interpolants^i"— , and the current- 
driven homogenization theory^li^. The latter approach 
has gained considerable traction lately2£"— . In this pa- 
per, we analyze this theory as an instructive case in point. 

One of the co-authors has already published^ a the- 
oretical analysis of the current-driven excitation model 
(not related to the theory of homogenization). However, 
since multiple claims have been made that the current- 
driven homogenization approach is rigorous, completely 
general and derived from first principles^ir— , it deserves 
additional scrutiny. Also, our previous analysis was 
mainly theoretical and no numerical examples were given. 
But the best test of any EMT is the test of its predictive 
power. It appears, therefore, useful to investigate the 
predictions of current-driven homogenization by using a 
simple exactly-solvable case of one-dimensional periodic 
medium. 

In fact, current-driven homogenization has been al- 
ready applied to such mediap^i^. However, the trans- 
mission and reflection coefficients T and i? of a layered 
slab have not been studied in these references. Instead, 
the nonlocal permittivity tensor S(a;,k) (defined below) 
was computed numerically. Current-driven homogeniza- 
tion of Refs. [24II25I entails an additional step in which 
S(a;,k) is used to compute purely local effective tensors 
e and /i (in non-ccntrocymmctric media, magneto-electric 
coupling parameters must also be introduced) and then T 
and R according to the standard formulas [e.g., see equa- 
tion (pij) below] . The nonlocal tensor E(w, k) can be used 
for this purpose only when complemented with additional 
boundary conditions (ABCs) whose nature and mathe- 
matical form are far from being clear; in any event, this 
computation has not been done. In addition, Refs. ISSllSa 
do not provide a closed- form expression for E(a;,k) not- 
ing that it is too involved^. 

In this contribution, we derive a closed-form expression 



for S(a;, k) in the case of s-polarization. Consideration of 
p-polarization is not mathematically difficult but is not 
needed for our purposes. We follow the current-driven 
homogenization methodology to derive closed-form ex- 
pressions for the local tensors e and fi. Then we use this 
result to compute T and R of layered slabs. In Sec.|ITl we 
summarize and discuss the prescription of current-driven 
homogenization of Refs. 24,25). In Sec. lIIII wc use this pre- 
scription to obtain closed-form expressions for the case 
of a one-dimensional layered medium. In Sec. II VI we list 
for reference the relevant formulas for the transmission 
and reflection coefficients of layered and homogeneous 
slabs. Numerical examples are given in Sec. [V] Here 
we compute the local effective parameters obtained by 
current-driven homogenization, by the S-parameter re- 
trieval method and by the classical (standard) homoge- 
nization approach. We then use these results to compute 
T and R and to compare the latter to the exact values for 
finite layered slabs. In Sec. IVIl we present a Bloch-wave 
analysis of current- driven homogenization. Secs. lVlIl and 
I Villi contain a discussion and a summary of the obtained 
results. Some technical details of the derivations and 
method used in this paper are given in the appendices. 



II. CURRENT-DRIVEN HOMOGENIZATION 

The current-driven homogenization theory is formu- 
lated for an infinite periodic medium and consists, essen- 
tially, of two steps. 

In the first step, one derives or computes numerically 
the nonlocal permittivity tensor E(w, k), which is defined 
as a coefficient between the appropriately averaged fields 
D(r) and E(r). The exact prescription for this computa- 
tion is given below. One could, potentially, stop at this 
point and attempt to use I](a;, k) directly to compute the 
physical quantities of interest. However, due to the ex- 
plicit dependence of S on k, it is not clear how to do so. 
At the very least, this computation would entail the use 
of ABCs. Since current-driven homogenization does not 
consider the physical boundary of a sample, derivation 
of the ABCs is outside of its theoretical framework. Be- 
sides, the use of the ABCs would defeat the very purpose 
of homogenization because all the applications of meta- 
materials discussed so far in the literature rely heavily on 
the existence of local constitutive parameters. 

Hence there exists a second step in which the non- 
local tensor I](cj, k) is used to derive purely local ten- 
sors e and fi (here we restrict attention to media with a 
center-symmetric lattice cells and do not introduce or dis- 
cuss magneto-electric coupling parameters). This second 
step is based on the proposition that, at high frequen- 
cies, magnetization of matter is physically and mathe- 
matically indistinguishable from weak nonlocality of the 
dielectric rcsponse^Sri^. We will give an exact prescrip- 
tion for completing this step, too. 

We now turn to the mathematical details needed to 
complete the two steps mentioned above. We work in 



the frequency domain and the time-dependence factor 
exp{—iujt) is suppressed. The dependence of various 
physical quantities on uj is assumed but not indicated 
exphcitly except in a few cases, such as in the notation 
S(w,k), where both arguments w and k are customarily 
included. The free-space wave number fco and wavelength 
Aq arc defined by 

fco = uj/c , Ao = 27r/fco . 

Finally, the Gaussian system of units is used throughout. 



A. Step One: calculation of the nonlocal 
permittivity tensor E(aj,k) 

Consider an infinite, periodic, intrinsically-nonmagne- 
tic composite characterized by the permittivity function 

etrue(r)- 

Here the subscript has been used to indicate that 
etrue(r) is the true spatially- varying parameter of the 
composite, as opposed to the spatially-uniform effective 
parameters e and fi. We assume for simplicity that the 
composite is orthorhombic so that 

etmc(a; + hx,y + hy, z + h^) ^ etruo(a;, y, z) , (1) 

where h^^hy and h^ are the lattice periods. Note that 
etruc(r) is a macroscopic quantity and that we consider 
the composite exclusively within the framework of macro- 
scopic electrodynamics. 

In the current- driven homogenization theory, it is as- 
sumed that the system is excited by an "impressed" or 
external electric current Jext (r) in the form of an infinite 
plane wave, viz. 



Jcxt(r) 



47ri 



-Je* 



(2) 



Here J is the amplitude, k is an arbitrary wave vec- 
tor which defines the "forced" Bloch-periodicity, and the 
Ld/Ani factor has been introduced for convenience. Note 
that Jcxt(r) is not subject to constitutive relations and 
is not equivalent to the current induced in the medium 
by the electric and magnetic fields. Maxwell's equations 
for the system just described have the following form: 

VxH(r) = -ifco[etruc(r)E(r)-fJe*'^-"-] , (3a) 

V X E(r) = ifcoH(r) . (3b) 

In some generalizations^^, a similar wave of magnetic cur- 
rent is included in pb[) . However, inclusion of electric 
current only will prove sufficient for our purposes. 

Obviously, the solution to ([3]) has the property of 
"forced" Bloch-periodicity^. This can be expressed 
mathematically as 



E(r) 



^ik-r 



F(r) 



(4) 



where F(r) satisfies the periodicity condition ([T]), and 
similarly for all other fields. The averaging procedure 



is then defined as "low-pass filtering" of the fields (e.g., 
Ref . [261) . The averaged quantities are defined according 
to 



- / e-*''E(r)dV = - / F(r)d3 

y Jc y Jc 



(5) 



Here V — h^hyhz — Jfjd^r and C denotes the unit 
cell. Similar definitions can be given for averages of all 
other fields, including the field of displacement D(r) = 
etruc(r)E(r). 

The nonlocal permittivity tensor is then defined as the 
linear coefficient between Dav and Eav, viz. 



Dav = S(^, k)Ea 



(6) 



If all Cartesian components of Eav and Dav arc known, 
^ contains three linear equations for the tensor elements 
of S(a;, k). Since the external current ([2]) is a convenient 
fiction rather than a mathematical representation of a 
physical object, we have complete freedom in choosing 
the direction of its amplitude J. By considering three 
different polarizations of J, we can construct a set of nine 
linear equations. However, in non-gyrotropic media, the 
tensor E(a;,k) is symmetrio^ and has, therefore, only 
six independent elements. We can force the set to be 
formally well-determined by requiring that k • J = 0. 

In this regard, it is useful to note that the aver- 
aged fields satisfy k-spacc Maxwell's equations with a 
spatially-uniform source^: 



k X Hav = -fc-0 (Dav + J) , k X Eav = ^'oH 



O-tlav 



(7) 



Consequently, k • (Dav -I- J) = 0. If k • J ^ [the current 
wave in ^ is not transverse], we also have k • Dav 7^ 0, 
which is a troubling occurrence. It means, of course, 
that, in addition to the external current ([2]), we have 
included into consideration an external wave of charge 
density Pcxt(r) = (k • J/w) exp(ik • r). However, in the 
homogenized sample, we expect V • D = to hold. In 
this paper, we use only a transverse external current wave 
but note that more general excitation schemes have been 
considered^^. 

Let us further specialize to the case of a two-component 
composite in which the function etruc(r) can take two 
discrete values ea and £&. We will write C ~ CaUCb and 
etiuc(r) = ea a r e Ca, etiuc(r) = £& if r e Cfc. In this 
case, Eav = Qa + Qfc, Dav = faQa + ^bQb, whcrc 

Q„ = / F(r)dV , Qb= f F(r)dV . 

■JCa JCb 

Therefore, equation ^ takes the form 

(QaCa + QbCfc) = S(W, k) (Qa + Qfc) • (8) 

From the linearity of ([3]), we have Qa = TaJ, Qfc = t^J, 
where Ta and rfc are two tensors. If Ta + Tfc is invcrtible, 
we can solve dSI to obtain 



S(a;, k) = {TaEa + net) {Ta + n) 



(9) 



The above equation implies that introduction of the ex- 
ternal current ^ is not really required to define the 
function S(w, k) mathematically. In fact, this statement 
is general and applies to any periodic structure in any 
number of dimensions, as long as the intrinsic constitu- 
tive laws are linear. In Sec. IVIi we will demonstrate the 
same point from Bloch-wave analysis. In Sec. IVII A| we 
will show that this function does not characterize the 
medium completely but can only be used to find the law 
of dispersion. 



B. Step T'wo: calculation of local parameters 

The proposition that magnetization (nontrivial mag- 
netic permeability) of matter is indistinguishable from 
nonlocality of the dielectric response is based on the 
equivalence of expressions for the induced current that 
are obtained in both models for infinite plane waves. Here 
we recount these arguments insomuch as they are needed 
for deriving the main results of this paper. 

Consider two electromagnetically-homogeneous media. 
The first medium is characterized by a nonlocal permit- 
tivity tensor S(a;,k) and fi = 1. In fact, the auxiliary 
field H is not introduced for this medium, so that fj, is, 
strictly speaking, not defined. The macroscopic electro- 
dynamics is then built using the fields E, B and D with 
the account of spatially-nonlocal relationship between D 
and E. The induced current in such a medium is given 

by 



J 



(1) 



ind 



47r 



[I](a;,k)-1]E 



(10) 



Here we assume, as is done in all relevant references, that 
E is an infinite plane wave with the wave vector k. 

The second medium is characterized by purely local 
tensors e and /i and the induced current in this medium 
is given by 






cV X M 



(11) 



where M is the vector of magnetization. Using the def- 
inition of M and macroscopic Maxwell's equations, we 
can also write 



J 



(2) 



ind 



47r 



(e-i) 



i"2 



kx (1 



^J■ 



kx 



E 



fl2) 



This expression can be compared to (|10p . In general, of 
course, there is no equivalence between pO)) and ([T2|. 
But in the so-called weak nonlocality regime [defined 
more precisely after Eq. ([T^ below], I](w,k) is well ap- 
proximated by its second-order expansion in powers of 
k. In this case, one can look for the condition under 
which P^ and ((T^ agree to second order in k. In non- 
gyrotropic media, the expansion of E(aj, k) has the form 



S(w,k) = S(a;,0) - -^k x /Jk x 



(13) 



where /3 is a tensor and the minus sign has been used for 
historical reasons. 

We are interested in the condition under which the 
expression in the square brackets in (J12p is equivalent 
to [S(w, k) — 1] computed to second order in k. It is 
easy to see that this condition is e = T,{lu,0) and fi = 
(1 — /3)~^ where the last equation implies tensor inverse. 
In isotropic media ^ and /3 arc reduced to scalars. In 
the case of cubic symmetry, when the tensors /3 and /i 
are diagonal in the rectangular reference frame XYZ, 
we have ^laa = (1 ^ l^aa)^^, where a = x,y, z. Note that 
all three principal values of /? can now be different. The 
conclusion that is typically drawn from this analysis^ is 
that the introduction of local parameter // is physically 
indistinguishable from the account of the second-order 
term in expansion (J13p . 

If the function S(a;,k) is known (it is computed di- 
rectly in Step One of the current-driven homogenization 
prescription), the tensor /3 can be easily computed from 
P3p . In the case of cubic symmetry, the relevant formu- 
las are 



Pxx 
Pyy ■ 



P /92y 



2 


dkl 


fi,g 


5^1],, 


2 


dkl 


ftp 


d^^xx 



dkl 



kg O ^zz 

Y~dkf 

kfj O L-izz 

T dkl 

P rP-V 
"-O '-' '^yy 

y dkl 



,2 " ^yz 

' °dkydkz 

2 d 2-ixz 

'°dhM'z 

2 i-l ^xy 



-k, 



(14a) 
(14b) 
(14c) 



All derivatives in the above equations must be evaluated 
at k = 0. 

We can now formulate the condition of weak nonlocal- 
ity more precisely. Let's assume that we have applied the 
prescription and computed the local parameters e and fi 
at a given frequency. These parameters can now be used 
to compute the natural wave vector of the medium, q, 
using the dispersion relation [e.g., for a uniaxial crys- 
tal, see p3p below]. We then evaluate the nonlocal per- 
mittivity I](a;,k) at k = q. The nonlocality is weak 
if the expansion ^3]) computed to second order accu- 
rately approximates E(aj,q). Thus, in the weak nonlo- 
cality regime, higher-order terms in the expansion (|13p 
can be neglected. 

At this point, we can make two important observa- 
tions. First, the above discussion applies only to infinite 
media. In any finite magnetic medium, additional surface 
currents exist. These currents are not included in (jlip . 
Consequently, the equivalence of currents is, in principle, 
not complete: it docs not apply to the surface currents. 
As a result, introduction of a nontrivial magnetic perme- 
ability and a dynamic correction to the permittivity^, as 
described above, can yield a first nonvanishing correction 
to the dispersion relation but not to the impedance of the 
medium. A related point is that Jjjjj(k) is not reduced 
to a quadratic form in finite samples even in the case of 
natural magnetics. 

The second observation is more subtle. The local pa- 
rameters that satisfy the requirement of current equiva- 




FIG. 1: Geometry of wave propagation in the case of s- 
polarization. Here k = {kx,0,kz) is the wave vector of the 
external current wave [Eq. ([2])]. A finite symmetric slab con- 
taining N = 6 unit cells is shown. Each cell consists of three 
layers of the widths {a/2, b, a/2). Equivalently, we can view 
the unit cells as consisting of two layers of the widths {a,b), 
provided that one half of the first a-type layer has been cut off 
and moved from the left face of the slab to its right face. Note 
that the sample shown in the figure has a center of symmetry. 



lence are not unique when the current is evaluated on- 
shell, that is, for k = q. There exists an infinite set of 
such parameters, related to each other by the transfor- 
mation (|37p (stated below), all of which yield exactly the 
same law of dispersion and the same induced current (|12p . 
However, in current-driven homogenization, the variable 
k in P^ is viewed as a free parameter. If we follow this 
ideology and require the current equivalence to hold for 
all values of k, we would obtain an unambiguous "pre- 
scription" for computing the local effective parameters. 
But the only physically-realizable case is k = q. There- 
fore, it is not clear why the pair of effective parameters 
predicted by current-driven homogenization is "better" 
than any other pair obtained by the transformation (1371) . 
This point will be illustrated numerically in Sec. IV Al 



III. EXACT SOLUTION IN 
ONE-DIMENSIONAL LAYERED MEDIUM 

The geometry considered is illustrated in Fig. [1] A 
one-dimensional periodic medium consists of alternating 
intrinsically nonmagnetic layers of widths a and b and 
scalar permittivities ea and eb, respectively. The period 
of the system is given hy h ^ a + b. The layered medium 
described here can be considered as a special case of the 
three-dimensional orthorhombic lattice obtained in the 
limit hx — hy ~ 0, h^ = h > 0. Note that the medium 
shown in Fig. [1] is finite and terminated by half- width 
a-type layers. However, in the current-driven homoge- 
nization theory, the medium is assumed to be infinite. 

We will consider only the special case of s-polarization, 
when the wave vector k lies in the plane XZ oi the rect- 
angular frame shown in Fig. [T] and the amplitude J of 
the external current ([2]) is coUinear with the F-axis, so 



that k — {kj:,0,kz) and J = (0, Jy,0). According to 
([T^. this is sufficient to uniquely define the following 
elements of the effective permittivity and permeability 
tensors: €xx ^ ^yyj l-^xx ^ l-^yy and ^zz- 

We will need to introduce the following notations: 



Kofca 


kb — k^eb , 




(15a) 


^kl 


~kl , Kb = -y/Zcg - 


- i-2 


(15b) 


Kad , 


(j)b = Kbb , 




(15c) 


kza , 


0b = kzb , 




(15d) 


a/h, 


Pb = b/h , 




(15c) 



Pa 



and also the standard homogenization result for a peri- 
odic layered medium: 



e|| ^Pa<^a +Pb^b 
Mil = Ai± = 1 • 



1 



Pa/<^a +Pb/^b 



(16a) 
(16b) 



Here the quantities indexed by "|j" and "_L" give the 
standard homogenization results for the elements of the 
permittivity and permeability tensors that correspond to 
the direction of the electric field parallel and perpendicu- 
lar to the layers, respectively. Throughout the paper, the 
branches of all square roots are defined by the condition 
TT < arg(yi) < 0. 

We wish to solve Eq. ^ in which etruo(r) is equal to 
Ca in the a-type layers and to Cb in the 5-type layers. 
Without loss of generality, we can consider the unit cell 
Q<z<h = a-\-b, which contains two layers: the first 
layer (a-type) is contained between the planes z = Q and 
z = a and the second (6-type) layer is contained between 
the planes z — a and z = h. We can seek the solution in 
each homogeneous region excluding its boundaries as a 
particular solution to the inhomogeneous equation plus 
the general solution to the homogeneous equation, viz, 

Ey{x, z) = e''^-- [gp{z) + Sg{z) ] , (17a) 

jy,(x,z) = e'''--[jrp(z) + ^4(z)] . (17b) 

Here the subscripts "p" and "g" denote the particular 
and the general solution, respectively, and the overall 
exponential factor exp(ik • r) is written out explicitly. 
The particular solution is given by 

4,(z) = jyf{z) , .ye^{z) = -^/yfi^) ' (18) 



where 



k"^ 



/(^) 



We emphasize that 
open intervals < z 



fa , < z < a 
= fb , a < z < h 

T8| is a particular solution in the 
< a and a < z < h. To satisfy 



/r2 
P - fc2 



boundary conditions at the interfaces, we must add to 
(fT5)) the general solution to the corresponding homoge- 
neous problem. The latter can be easily stated: 



A^) = 



JyAe-^'^^'Fsiz) , 



J^Az) = -j^^JyAe-^'^'^FMz) 



(19a) 
(19b) 



where 



and 



Fsiz) 



A = fb-fa = 



ftp 






1,2 _ 1,2 jL2 _ jL2 



(20a) 



A^e^"-^ + Bac-^"-'' , < z < a 
e'''"[-A6e*'^''(^-'^) - She-^'^o^^-'^)] , a<z<h ' 

F^(z) = (20b) 

— e*^-[-Ae*''''(^-'') + Bte-''"'^^-''^] , a<z <h 

In these expressions, various z-independent factors have 
been introduced for convenience and Aa, Ba, At, Bt is a 
set of coefficients to be determined from the boundary 
conditions. The latter require continuity of all tangential 
field components at the interfaces z = 0,a,h and can be 
stated as follows: 



F^(a-0)-F^(a + 0) = e'''» 
i^^(a-0)-F^(a + 0) = e" 



JBa 



JBa 



F^iO) - 
FM^) 



'Fg{h) = 
"'F^ih) 



JBa 



(21a) 
(21b) 
(21c) 
(21d) 



where 

^ = cos{kzh) - cosiq^h) (23) 

and the closed- form expressions for £/a, ^a, ■^b and ^b 
are given in Appendix [X] In (|23p . Qz is the z-projection 
of the natural Bloch wave vector q computed under the 
assumption that its X-projection is equal to k^ (that is, 
Qx = kx). The factor cos{qzh) is defined by the equation 



cos{qzh) 



-- COs{(t)a) COs{(j)b) 



sm 



,)sin(0,). (24) 



Evidently, if kz = ig^ + 2Tm/h, where n is an arbitrary 
integer, the matrix in ()2ip is singular. 

We now simplify the expression for the electric field 
Ey{x, z). After some rearrangement, we can write 

Ey{x,z) = 
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The yy-componcnt of the nonlocal permittivity tensor 
E(a;,k) is computed by using (l8|) or (j9|), which results in 



Eyy(w,k) 



Qafa + Qb^b 



(26) 



This results in a set of four equations for the unknown co- 
efficients Aa, Ba, Ab, Bf,, which are stated in Appendix El 

It may seem confusing that the right-hand side in ((2T|) 
does not go to zero when Jy — > 0; in fact, (PT|) does not 
contain Jy at all. However, the electromagnetic fields 
Ex{x,z) and Hy{x,z) computed according to (|17p - (fT9|) 
are proportional to Jy. Note that the most general solu- 
tion to ([2]) is a superposition of the solution derived here 
(whose amplitude is proportional to Jy) and the natural 
Bloch mode of the medium with an arbitrary amplitude. 
To remove the nonuniqueness, one can either consider the 
boundary of the medium and thus abandon the infinite 
medium model, or, alternatively, apply the additional 
boundary condition requiring "forced" Bloch-periodicity 
(|^. The latter approach is used in current-driven ho- 
mogcnization and in the derivations of this section. 

The solution to (|2ip is given by 



where 



A„ 



2& 



Bn 



2<3i 



Au 



221 



Bh 



^b 

22 



(22) 



Qa -= / Faiz)dz , Qb= Fb{z)dz 



(27) 



The integrals in (|27|) arc easily computed analytically; 
these intermediate results are omitted. Note that Q^ and 
Qb depend implicitly on both lo and k, which are con- 
sidered as mathematically-independent variables in the 
current-driven homogenization theory; this dependence 
is indicated explicitly in the notation E(cj, k). 

Equations (P5|) - (P7)) together with the expressions for 
the expansion coefficients given in Appendix |X] consti- 
tute a closed- form solution for I]j,j,(w,k). This solution 
contains only elementary functions, has no branch am- 
biguities [see the note after Eq. ([T5|)] and can be easily 
programmed. Note that the quantities Qa and Qf, defined 
in (|27p have no singularities when viewed as functions of 
k. However, I]j,j,(aj,k) has singularities at the roots of 
the equation Qa + Qb = 0. This completes Step One of 



the current-driven homogenization prescription, at least 
for the case of s-polarization. 

We now proceed with Step Two. For the local effective 
permittivity, we have 



'^vv 



^yyiuj,0) 



Qa^a + Qb^b 



Qa + Qb 



k=0 



We note that this expression contains the dynamic cor- 
rection to the permittivity^. From symmetry, we also 



have 



^yy 



The remaining nontrivial component 



of the permittivity tensor is e^^; this element cannot be 
computed by considering only s-polarization of the exter- 
nal current. 

To compute the elements of the permeability tensor, 
we use the first equality in (|14ap and the second equality 
in (|14cp . More specifically, we have 



^2 fl2y 



k=0 



' '^"^ 2 dkl 



k=0 



Using (pS)) . we can obtain the following formulas for P^x 
and j3zz in terms of Qa and Qf,: 



Jzz) 



Pxx — f^o\^b ^a) 



^a^h ^b^a 



(zz) 



{Qa+QbY 



Pzz = fco(£b ^ ^a) 



QaQt^ 



QbQ^r'' 



[Qa + Qb? 



k=0 



k=0 



(28a) 
(28b) 



Jxx) 



In these expressions, Qa denotes the second derivative 
of Qa with respect to kx evaluated at k = 0, etc. In 
deriving (pS)) . we have used the fact that Qa' = Qa' = 

Qb ~ Qb ~ ^- Note that equation (^5]) is invariant 
under the permutation of indexes a -^r^ b. 

The elements of the effective permeability tensor are 
expressed in terms of Pxxi Pzz as 



1 



1 



l^XX 



l-Pz 



(29) 



From symmetry, we also have fXyy = jjlxx- Closed-form 
expressions for exx = ^yy, ^^xx = f^yy and /Jzz are given 
in Appendix |B] These expressions contain only elemen- 
tary trigonometric functions but are fairly cumbersome. 
However, the small- h asymptotic approximations of these 
expressions have the following simple form: 



- ^yy - H + -^ 

X {PaPbf {kohf + 0{h^) 



l^xx 



^iyy = 1 + 



(ea - Efo) 



240 

X (PaPbf il+'2PaPb)ikohy 



om 



(ea -efc)^ 



720 



X {paPb) (1 + 2paPbKkohf + 0{h^) 



(30a) 



(30b) 



(30c) 



Thus, the first nonvanishing corrections to the effective 
permeability tensor are obtained to fourth order in h. 
Moreover, the corrections to fixx and ^izz differ by the 
constant factor —3. Consequently, current-driven homog- 
enization, when applied to the ID periodic structure con- 
sidered in this section, guarantees that at least one of the 
principal values of the permeability tensor has a negative 
imaginary part for sufficiently small values of h, provided 
that Im(ea - e^f ^ 0. 



IV. TRANSMISSION AND REFLECTION BY 
HOMOGENEOUS AND LAYERED SLABS 

In what follows, we will need to refer to the formulas 
for the transmission (T) and reflection [R) coefficients 
of homogeneous and layered slabs. These formulas are 
well known and are adduced here mainly for reference. 
However, they also reveal some important features that 
will help us analyze the numerical results of the next 
section. We still work in the geometry of Fig. [T] and, in 
analogy to (jlSbl) . denote the z-projection of the incident 
wave vector by kq, so that 



Ko = \Jkl - kl . 

Note that, for kx > fco, the incident wave is evanescent. 
All formulas given below arc parameterized by kx- 

Any slab with a one-dimensional distribution of elec- 
tromagnetic parameters is completely characterized by 
its transfer matrix M. Suppose the slab occupies the re- 
gion < z < L. Then the tangential components of the 
electric and magnetic field at the left and right faces of 
the slab are related by 



El 
Hl 



Mil Mi2 
M21 M22 



Eo 
Ho 



The general property of all transfer matrices is unitarity, 
from which it follows that det(Af) = 1. If the slab has a 
center of symmetry (as is the case in this paper), the left 
and right incidence directions are equivalent which can 
be mathematically stated as Mu = M22- Under these 
conditions, the transfer matrix can be written as^ 



M 



cos 61 (-i/^)sine' 
-i^sinO cos 9 



where 9 and 2f are the optical depth and the generalized 
impedance of the slab^. 

The transmission and reflection coefficients can be ex- 
pressed in terms of 9 and 3f as 



T - 
R 



1 



cos9~iX+{%,3f)siii9 
-~iX.{3fo,3f)sm9 



cos 9 — iX^ 



Wo, ^) suit 



(31a) 
(31b) 



where 



«3;,^.).i(|±f 



8 



and ^ = Ko/fco is the generalized impedance of free 
space (we assume that the slab is embedded in a vacuum 
or air). The quantities T and R defined in (|31]) relate 
the amplitudes of the transmitted and reflected tangen- 
tial field (electric in the case of s-polarization or magnetic 
field in the case of p-polarization) measured at the planes 
z = L (for T) or z = (for R) to the amplitude of the 
incident wave at ^ = 0. The specific expressions for 9 
and ^ depend on polarization and it should be kept in 
mind that, at normal incidence, Tg = Tp but Rg = —Rp, 
where the subscripts indicate the particular mathemati- 
cal expression applicable to a given polarization state. 

Specific expressions for 9 and ^ for homogeneous and 
layered slabs are given below. 

a. Homogeneous Anisotropic Slab. Consider a slab 
characterized by purely local diagonal tensors e = 
diag(e^,e_L,e||) and ^ = diag(/i_L, /i_L, ^n). Then 



where 



QzL , ^ = q;:/kar 



<lz = ^k^HI-Hl -kx{v\\/vi-) 



(32) 



(33) 



and 7] refers to fi for s-polarization and to e for p- 
polarization. 

b. Layered Slab. Consider a layered slab of total 
width L ~ Nh containing N unit cells arranged as shown 
in Fig. [T] Each cell consists of three consecutive layers 
of the widths {a/2, b, a/2), where a + b = h, and the 
permittivities ea and et, respectively. Then 



9 = q,L = Nq,h , 

sin (f>a cos (f>h — X_ sin 0h + X-|_ cos (f>a cos 0(, 



sin (f>a cos (/);, + X_ sin ( 



X^ cos (/)a COS ( 



(34a) 



(34b) 



Here <j)a,<Pb are defined in (|15cp . X± ~ X±{2fa, 2fb), 3fa 
and 3(b are the generalized impedances of each layer and 
Qz is the natural Bloch wave number of the medium de- 
fined by the following equation: 

cos{qzh) = cos{(j)a) cos((/)(,) - X+ sm{(j)a) sin(0h) . (35) 

Note that ([24|) is a special case of ([35]) (for s-polarization 
and nonmagnetic layers). Also Eq. ([35]) defines cos{qzh) 
but not sin(gz/i). The latter quantity can be computed 
by using one of the formulas 



sm{qzh) 



6^ 

— (sin (^a cos 0b + X_ sin (/)^, + X+ cos (^a 



sine 



= -^ (sin (i)a cos 0t - X_ sin 0,, + X+ cos 0^ sin 0b) , 

where SH" is determined by taking an arbitrary branch of 
the square root of (|34bp ; the resultant transfer matrix is 
invariant with respect to this choice. 



The formulas given above illustrate several impor- 
tant points. First, the transmission of a thick, highly 
transparent slab is very sensitive to small errors in 
qz. This is because the trigonometric functions such as 
cos 9 = cos{qzL) incur a substantial phase shift when qz 
is changed by ^ n/L. If L ^ cx) (and losses can still be 
ignored), any homogenization theory is expected to be 
unstable numerically because a small error in determin- 
ing the EMPs can propagate to become a significant error 
in T. This instability, however, is of little practical im- 
portance because, in most cases, the illumination is not 
monochromatic. What we are discussing here are, essen- 
tially, the resonances of a Fabry-Perot etalon. In most 
applications to optical imaging and microscopy, illumi- 
nation is more broadband than a line of a single-mode 
laser and the interference effects are unobservable. Under 
these conditions, the expressions pip can be regularized 
by Gaussian integration with respect to fco or L. How- 
ever, in this paper, we only consider strictly monochro- 
matic light. 

Second, an error in the generalized impedance will also 
result in an error in both T and R and, in some cases, 
this error can also be dramatic. An illustrative example 
is the case X+ — > — 1, which is the operation regime of 
the so-called perfect lens^. Indeed, we can re- write (|31ap 
as 



T = 



(1 - X+) cxp{i9) + (1 + X+) exploit 



(36) 



li X+ is exactly equal to —1, ([55]) predicts T = exp{—i9). 
For evanescent waves and a macroscopically-thick slab, 
the factor exp(i9) is exponentially small. Therefore, if 
we make a small error— in X+, say, X^ = —1 + 6, 
such that \exp{i9)\ < \S\ < 1, Eq. dMl) would predict 
T ~ {2/S) cxp{i9), which is dramatically different from 
the former result. We note that there are other similar 
situations and that the condition Xj^ w — 1 is not special 
in this respect. 

Therefore, it is incorrect to assume that the main goal 
of a homogenization theory is to predict correctly the op- 
tical depth: both the optical depth and the impedance 
are important. But can local EMPs be found in such a 
way as to predict 9 and J^ correctly and simultaneously? 
Standard EMTs do allow this asymptotically in the limit 
/i — > 0. In contrast, extended EMTs that consider an in- 
finite medium cannot predict correctly the impedance of 
the medium. The current-driven homogenization theory 
is of this variety: it predicts correctly the first nonvan- 
ishing correction to 9 (compared to the standard homog- 
enization result) but does not provide any meaningful 
corrections or approximations to 3f. Moreover, this cor- 
rection to 9 yields a valid approximation only in a limited 
range of ft,, as will be shown below. In a more general 
case, it is not possible to find local EMPs that predict 
correctly (for all angles of incidence) even 9 alone. There- 
fore, the claims that the current-driven homogenization 



theory is rigorous and completely genera' 
gerated. 



24-26 



are exag- 



We finally note that the above analysis could be ex- 
tended to three-dimensional orthorhombic lattices by in- 
tegrating out higher-order harmonics in the ccy-plane, i.e. 
by low-pass filtering. 



V. NUMERICAL EXAMPLES 

In this section, we consider several examples of com- 
puting the local effective parameters of one-dimensional 
layered media according to the current-driven homogc- 
nization prescription. We note that the maximum pos- 
sible value of the numerical factor (paPb) (1 + '^PaPb) in 
(PH]) is 3/32 and it is achieved when p^ = pb = 1/2. These 
volume fractions are used in all the numerical examples 
shown below. The 6-type medium is assumed to be air 
or vacuum with eb = 1. For the a- type medium, three 
different examples will be considered. In Example A, the 
a-type medium is a lossy dielectric considered at a fixed 
frequency and varying values of h and kx- Example B is 
similar to example A, but the a-type medium is a con- 
ductor. In Example C, the a-type medium is an idealized 
Drudean metal considered at a fixed value of h, k^ = 
and varying Aq. Thus, in Example C we account for the 
frequency dispersion in metal. 

The results of current-driven homogenization will be 
compared to the standard homogenization result ()16|) and 
to the results of S-parameter retrieval^>^"— . As is well 
known, retrieving effective parameters of a slab from the 
transmission and reflection coefficients at normal or near- 
normal incidence is an ill-posed inverse problem. We 
have used several different methods of regularizing the 
inverse solution, some of which have been proposed in 
the literature^ and others have been devised by us; see 
Appendix [C] for full details. All these various modalities 
of the retrieval technique yield approximately the same 
result when h/Xo <S^ 1 but deviate strongly for larger val- 
ues of this ratio. In the figures below, the retrieval results 
are shown only for the range of h/Xo within which the 
technique is numerically stable. 



A. Example A 

We start with the case where the a-type medium is 
a lossy dielectric characterized by Ca = 4.0 -I- 0.1 i at a 
given wavelength Aq. In Example A, we assume that h 
and kx can vary while Aq is fixed. Then the physically- 
measurable quantities of interest are functions of the di- 
mensionless variables h/Xa and kx/ko, and the actual 
value of Ao is unimportant. The sample consists of 
A'' = 50 symmetric unit cells of the type (a/2, &, a/2) 
arranged as shown in Fig. [TJ As noted above, we take 
a = b = h/2 in all numerical experiments. 

In Figs. [2]|4l we illustrate the predictions of current- 
driven homogenization for all the components of the per- 
mittivity and permeability tensors that can be obtained 
in s-polarization. The results are compared to the predic- 



tions of the S-paramcter retrieval method. It can be seen 
that current-driven homogenization produces the stan- 
dard homogenization result when /i — > 0. This much 
could be inferred from considering the asymptotic expan- 
sions (pO)) . In fact, for h/X^ < 0.2, all curves displayed 
in the figures are close to the standard homogenization 
result. We note that this is, approximately, the same 
range of h in which S-paramctcr retrieval is numerically 
stable. The interpretation of this fact is obvious: for 
sufficiently small ratios of h/Xo, the transmission and re- 
fiection properties of the sample are well fitted by purely 
local effective permittivity e and /x = 1. 

Both current-driven homogenization and S-paramcter 
retrieval provide corrections to the standard homogeniza- 
tion result. As long as these corrections are small, they 
can yield a homogenization result that appears to be 
"reasonable". Yet, all the dramatic features of current- 
driven homogenization occur for h/Xo > 0.2, where the 
S-parameter retrieval is unstable. In particular, the res- 
onance in fixx occurs at h/Xo ~ 0.4. Is this result phys- 
ically reasonable? We will address this question now by 
considering the transmission and reflection coefficients of 
a finite slab. 

In Fig.[5l we plot jTp and |_Rp at normal incidence as 
functions of h/Xo- It can be seen that the different meth- 
ods used to compute the transmittance and reflectance 
yield very similar result for h/Xo ^ 0.2 but very differ- 
ent results for h/Xo > 0.2. Overall, when both T and R 
are considered, current-driven homogenization does not 
provide a meaningful correction to the standard homoge- 
nization result P^ . In other words, at small h/Xo, both 
methods predict approximately the same result while at 
larger values of h/Xo both methods simultaneously fail. 
This is clearly visible in the case of R but is also true for 
T, which is very small when h/Xo > 0.4. In the latter 
case, both current-driven and standard homogenization 
generate relative errors in jTp of many orders of magni- 
tude, as could be veriflcd by utilizing logarithmic vertical 
scale (data not shown). Of course, this result is expected 
for standard homogenization, which is an asymptotic the- 
ory. However, the current-driven homogenization theory 
was claimed to have predictive power beyond the limit of 
small h and, in particular, in the region of the parameter 
space where it predicts nontrivial magnetic effects. This 
claim appears not to be supported by the data of Fig. \5\ 

Nevertheless, if we focus on T alone, current-driven 
homogenization provides a slightly more accurate result 
compared to standard homogenization when h/Xo is in 
a small vicinity of 0.2. Let us, therefore, consider in 
more detail transmission and reflection by the slab at 
h/Xo = 0.2. The EMPs obtained at this value of h/Xo by 
current-driven homogenization are listed in Table U and 
the dependence of |Tp and |i?p on the angle of incidence 
is illustrated in Fig. [51 In the case of T, current-driven 
homogenization provides a noticeable improvement over 
the standard homogenization result when kx < fco (in 
fact, the standard formula predicts the phase of T incor- 
rectly in this range of kx) but not when kx > ko- In the 
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FIG. 2: Example A. Real (top) and imaginary (bottom) parts 
of eyy as functions of h/Xo. The various curves shown are 
obtained as follows: CD - by current-driven homogenization 
[formulas given in Appendix[B]; AS - the small-/i asymptotes 
of the former [Eq. (|30p ]: RET - by S-parameter retrieval [see 
Appendix [C]; and AV - the standard homogenization result 
ey [defined in Eq. (|16ap ]. Insets show the details of all curves 
for small h. 



case of |i?|^, no improvement is observed. 

We will discuss below the reason why current-driven 
homogenization predicts jTj^ more accurately than the 
standard homogenization result at /i/Aq = 0.2, but it 
is useful to note right away that this has nothing to do 
with an accurate prediction of /.j. In fact, the values of 
fi computed by current-driven homogenization are nei- 
ther physically-significant nor optimal. To illustrate this 
point, consider the data of Fig. [7l Here we plot the real 
parts of T and R and introduce two additional curves. 
The first of these curves (labeled CD- REN) was obtained 
by taking the current-driven EMPs listed in Table U and 
rcnormalizing them according to the formula 



FIG. 3: Example A. Same as in Fig. [2] but for fixx- 
standard homogenization result Hxx = 1 is not shown. 



The 



e ^- ^e , n^ /i/$ 



(37) 



with the renormalization factor ^ = fi^x- In (|37p . e nd 
fi are tensors while ^ is a scalar. Renormalization (j37p 
does not affect the equivalence of the induced currents 
pni and (|12p when each formula is evaluated on-shcll, 
that is. for k = q. The renormalizcd parameters are also 



given in Table HI It can be seen that the renormalizcd 
EMPs have dramatically different values of /x — 1. Yet, 
T and R computed by using both sets of parameters are 
virtually indistinguishable. 

Moreover, the EMPs obtained by current-driven ho- 
mogenization are in no way optimal if the goal of homog- 
enization is to fit the transmission and refiection data as 
closely as possible. The latter aim is, in fact, achieved by 
the S-parameter retrieval procedure. In Fig. [7] we show 
an additional curve (labeled RET), which was computed 
using EMPs obtained by S-parameter retrieval. More 
specifically, eyy and ^^x have been computed by Method 
2 and ^^z was computed by Method 3, where the various 
methods of S-paramcter retrieval are described in Ap- 
pendix [Cl The particular choice of methods is explained 
as follows: Method 2 is more stable numerically but, un- 
like Method 3, it does not allow one to compute Hzz- 
Returning to Fig. [71 we observe that the curve labeled 
RET provides a much better fit to both T and i? in a 
wide range of incidence angles. This is in spite of the 
fact that the EMPs labeled as RET in Table |T] are very 
different from those labeled as cither CD or CD-REN. 

At this point, wc note that the magnetic effects pre- 
dicted by current-driven homogenization at /i/Aq = 0.2 
arc tiny; |/i — 1| does not exceed ^ 0.01 and the condi- 
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FIG. 4: Example A. Same as in Fig. [2] but for /i^z- The 
standard homogenization result n^z = 1 is not shown. 



TABLE I: EMPs for Example A obtained by various meth- 
ods at h/Xo = 0.2 and h/\o = 0.3. AV - the standard 
homogenization result; CD - current-driven homogenization; 
CD-REN - current-driven homogenization with renormaliza- 
tion (|37|l : RET - retrieved parameters. The numbers have 
been rounded off to three significant figures. However, all 
plots shown in this paper utilize either double or quadruple- 
precision computations. 



tion of weak nonlocality is very well satisfied. Yet the 
relative errors in T and R produced by current-driven 
homogenization are significant - they are at least of the 
same order of magnitude as |// — 1| or greater. In par- 
ticular, the relative errors in Ke{R) or \R\'^ at normal 
incidence exceed 100%. To distinguish the two effects, 
it would suffice to measure the reffection coefficient at 
normal incidence. 

Now let us turn to the case H/Xq — 0.3, which is il- 
lustrated in Figs. I8I9I S-parameter retrieval is unstable 
at this point and standard homogenization is inapplica- 
ble (the corresponding curves are not shown). Therefore, 
if current-driven homogenization could produce reason- 
able predictions at h/Xo = 0.3, it would constitute a 
valid and useful approximation. However, the data of 
Figs. I8I9I do not support this hypothesis. The relative 
errors in R and T for current-driven homogenization arc 
at this point dramatic and can be many orders of mag- 
nitude. For kx < fco, the relative errors are many orders 
of magnitude in T and at least an order of magnitude in 
R. Both the phase and the amplitude of T and R are 
predicted with significant errors in the whole range of k^ 
considered. We note that, in the bottom plot of Fig. [HI 
the real part of R is reproduced correctly at normal in- 
cidence by current-driven homogenization. However, as 



can be seen from the data of Fig. [51 |i?p at normal inci- 
dence is predicted with a large error. Consequently, Imi? 
is also predicted with a large error (data not shown). 

All this is in spite of the fact that current-driven ho- 
mogenization does not yet predict any dramatic magnetic 
effects at h/Xo = 0.3. Indeed, for this value of h/Xo, 
current-driven homogenization predicts that |/^— 1| does 
not exceed « 0.1. Moreover, in can be seen from the data 
of Figs. 1819] that renormalization of the EMPs according 
to ((37)) does not worsen or, indeed, noticeably modify the 
predictions of current-driven homogenization. 

At even larger values of h/Xo, predictions of current- 
driven homogenization are widely inaccurate. In par- 
ticular, current-driven homogenization cannot be relied 
upon at h/Xo ~ 0.4. when fi^x experiences a dramatic 
resonance. This is evident from the data of Fig. [S] and 
there is no need to support this conclusion with addi- 
tional graphics. A question, however, remains: why did 
we observe a moderate improvement in T at h/Xo = 0.2? 
The reason is that current-driven homogenization pro- 
vides a first nonvanishing correction to the Bloch wave 
number q^ but not to the impedance 3f. Both quantities 
enter the formulas for T and R (|5T|) . As was discussed in 
Sec. IIVI under some circumstances, T can be more sen- 
sitive to errors in q^ than to errors in i2°. This does not 
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FIG. 5: Example A. Absolute values squared of the trans- 
mission (top) and reflection (bottom) coefficients at normal 
incidence as functions of h/\o. The various curves shown are 
obtained as follows: EX - the exact result [Eqs. pip ]: CD - 
equivalent homogeneous slab with current-driven EMPs; and 
AV - same as above but for the EMPs obtained by standard 
homogenization. 



mean that errors in i2° are insignificant. As could be seen 
in Figs. [3 an error in the impedance that current-driven 
homogenization entails translates into an error in T and 
R, which is at least of the same order of magnitude or 
greater than |^ — 1|. 

The above point is illustrated in Figs. IIOIIII Here we 
plot q^h and i2° at normal incidence as functions of h/Xg; 
predictions of current-driven and standard homogeniza- 
tion are compared to the exact result given in ([M|) . (|35p . 
It can be seen that, at /i/Aq ~ 0.2, current-driven homog- 
enization provides a slightly more accurate result for q^. 
In the case of 3f, current-driven homogenization does not 
provides a better approximation at any value of H/Xq. 



B. Example B 

We now turn to the case when the a-type medium is a 
high-conductivity metal with ea = —3 + O.Oli at a given 
fixed wavelenghth Aq. The sample consists of 5 symmet- 
ric unit cells of the type (a/2,6, a/2), where a = 6 as 



FIG. 6: Example A. Absolute values squared of the trans- 
mission (top) and reflection (bottom) coefficients computed 
as functions of k^/ko for h/Xo = 0.2. Same curve labels as in 
Fig. [5] have been used. 



before. EMPs for Example B are plotted as functions 
of h/Xo in Figs. [T^lfHl Current-driven homogenization 
predicts that n^x experiences a resonance near the point 
h/Xo = 0.75 while fizz exhibits no dramatic effects. 

In Fig. [151 we display the predictions of various the- 
ories for |Tp and |i?p. The conclusion that can be 
made is that current-driven homogenization does not pro- 
vide a meaningful correction or a noticeable improve- 
ment of precision compared to standard homogenization 
in the whole range of H/Xq considered. In fact, there 
are fairly significant intervals of h/Xo (clearly visible in 
the insets) in which standard homogenization predicts 
correctly |Tp « and |i?p ~ 1 while current driven ho- 
mogenization is widely off the mark. Perhaps, current- 
driven homogenization can be credited with predicting 
a transparency window which exists in reality for rela- 
tively large values of h/Xo and which is not predicted for 
obvious reasons by standard homogenization formulas. 
Unfortunately, the transparency window is predicted for 
wrong values of h/Xo and, in the true transparency win- 
dow, both T and R are predicted with the wrong phase 
and amplitude. This is illustrated in Figs. (TH] and [T7] 
where we plot real and imaginary parts of both T and R 
as functions of k^/ko. It can be seen that there is no cor- 
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FIG. 7: Example A. Same as in Fig. |6]but for the real parts 
of T and R. Additional curve labels; CD-REN - current- 
driven homogenization with renormalization (JSTJI; and RET 
- S-parameter retrieval. See Table |T] for numerical values of 
the EMPs used for each curve. Note that CD and CD-REN 
curves are visually indistinguishable; EX and RET curves are 
indistinguishable in the upper plot but slightly different in the 
bottom plot. 



respondence between current-driven homogenization and 
exact result. 

The discrepancy is even more pronounced for the val- 
ues of /i/Ao such that the exact transmission coefficient is 
close to zero but current-driven homogenization predicts 
significant transmission. 



C. Example C 

In this example, vife consider spectral dependencies of 
|i?P and |rp at normal incidence vifith the account of 
frequency dispersion in the constituents of the compos- 
ite. The a-type medium is an idealized Drudean metal 
described by the permittivity function 



FIG. 8: Example A. Same as in Fig. E] but for /i/Aq = 0.3. 



eo 



u(uj + 17) 



(38) 



and the 6-typc medium is vacuum or air. The sample 
consists of iV = 10 symmetric unit cells of the type 



(a/2,6, a/2), where, again, a = b. We have chosen the 
parameters in ([55]) to represent the experimental values 
for silver: eo = 5 and Wp/7 = 500. The lattice period h is 
fixed so that h/Xp — 0.2, where Xp = lixLOpj c is the wave- 
length at the plasma frequency lOp. In the case of silver, 
Xp = 136nm so that h k, 27nm. The free-space wave- 
length Aq is varied. In this case, all physical quantities of 
interest can be expressed as functions of the dimension- 



less variables /i/Aq, 
in this example). 



Ao/Ap and k^jh^ (we will take k,^ 







We do not display the EMPs obtained by different 
methods as nothing qualitatively new compared to the 
previously considered examples emerges in Example C. 
Note that the current-driven permeability [i^x experi- 
ences a sharp resonance at /i/Aq ~ 0.38 while e^y and 
[izz do not exhibit any dramatic features. In Fig. I18[ 
we plot |rp and |i?p for Ap/Ao varying from to 2.5. 
The corresponding parameter /i/Aq varies from to 0.5. 
Again, the data clearly demonstrate that current-driven 
homogenization does not provide a meaningful correction 
to the standard result. 
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FIG. 9: Example A. Same as in Fig. [7] but for h/Xg = 0.3. 



VI. BLOCK- WAVE ANALYSIS OF THE 

CURRENT-DRIVEN HOMOGENIZATION 

THEORY 

In this section, we consider the current-driven homog- 
enization theory from a more general point of view. We 
assume that the medium is intrinsically-nonmagnetic, 
three-dimensional and periodic, and that its permittivity 
function etruc(r) satisfies the periodicity condition ([T|). 
Any such function can be expanded into a Fourier series 



et: 



ruc(,rj — / ^ EgG , 



(39) 



FIG. 10: Example A. Real (top) and imaginary (bottom) 
parts of the unit cell optical depth parameter, q^/i, at normal 
incidence as a function of h/Xo- 



pop , wc can find the relation between the expansion co- 
efficients Dg and Eg: 



Dg = ^ epEg_p . 



(40) 



Upon substitution of the expansions into ([3]) , we find the 
following system of equations for Eg: 



where 



, xn,r yn,, zti, 
g = 27r— -^ + :^ + — ^ 



are the reciprocal lattice vectors, which can be viewed as 
three-dimensional summation indices, and n^, riy and riz 
are arbitrary integers. We can seek the solution to ([3]) in 
the form of a Bloch wave 

E(r) = ^Ege*('^+s)-"- . 
g 

The displacement D(r) = etruc(r)E(r) can be similarly 
expanded. Given the periodicity of etruc(r) expressed in 



(k -f g) X (k + g) X Eg 



2^epEg_p + J(5go 



(41) 



Equations of this kind are well known in the theory of 
photonic crystal a^^'^^ , except that here we have included 
the free term JJgo. However, the following analysis (pro- 
posed by us earlier—) is rarely used. 

Let us write ((4T|) for the special cases g = and g 7^ 
separately. We note that Eg = Eav, where the low-pass 
filtered averages are defined in ([5]), and eo is the usual 
arithmetic average of etruc (without low-pass filtering). 
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FIG. 11: Example A. Real (top) and imaginary (bottom) FIG. 12: Example B. Real (top) and imaginary (bottom) 
parts of the generalized impedance 2f, at normal incidence as parts of eyy as functions of h/Xo- Same curve labels as in 
a function of h/Xo- Fig. [2]liave been used. 



Wc thus obtain: 
g = : 



k X k X Eav + fcg 



p#0 



0, 



gT^O : 



(42a) 



(42b) 



(k 4- g) X (k + g) X Eg 



Kn 



E 

PT^g 



fpt'g-p 



^gJ-Jav 



Note that only the first of these two equations is affected 
by our choice to include the external current in ([3|). We 
now utilize the linearity of equation (j42b|) . from which 
we can write: 



Pt^O 



epE_p = [S(w,k)-eo]Ea 



(43) 



where S(a;,k) is a tensor to be determined by solving 
(|42b[) . Here the factor eo (the average permittivity of 
the composite) has been introduced for convenience and 
does not result in any loss of generality. Then the g = 
equation (|42a|) takes the following form: 



k X k X Eav + fcg P(a;, k)Eav + J] = 



(44) 



This is, essentially, the same equation as ^. Conse- 
quently, I](aj, k) is the same tensor as the one appearing 
in the current-driven homogenization theory. 

We note that inclusion into Maxwell's equations of the 
external current ([2|) is not needed to compute I](aj,k), 
which is completely defined by the infinite set of equa- 
tions (|42bp . We can refer to this set as to the cell prob- 
lem. In what follows, we assume that the cell problem 
can be solved by means of linear algebra and that the 
tensor I](a;,k) can be computed. 

Since I](a;,k) is defined completely by solving the cell 
problem, it is useful to consider what would happen if 
we set J = in (|44l) . Obviously, this would result in an 
eigenproblem 



[(kx kx) + k^j:{uj,k)] Eav = 



(45) 



The above equation has nontrivial solutions only when 
k = q, where the Bloch wave vector q is determined 
from the equation 



det [(k X kx ) + fc^I](w, k)] = 



(46) 



The solution to this equation, viewed as a function of 
frequency, yields the dispersion equation of the medium. 



16 





Re/ijjj; 














1 








~ 






80 


— 




~-.^ 










0.9 


_ 




\ 






40 






1 



















40 





0.1 


0.2 1 




CD 

AS 








r 




RET 








SO 


1 


1 


1 1 1 




Ira/ixx 


1 


1 1 

1 


40 


_ 










0.0002 
















.'- 


' 









— 


_.--■' 


"^^ 






20 






-0.0002 




1 
























0.1 


0.2 




1 1 1 1 /^/^o 




0.2 



0.4 



0.6 



0.8 



0.2 



0.4 



0.6 



0.8 



FIG. 13: Example B. Same as in Fig. [12] but for /i^^. The 
standard homogenization result [jl^x = 1 is not shown. 



FIG. 14: Example B. Same as in Fig. [Tg] but for ^2^. The 
standard homogenization result /Xzz = 1 is not shown. 



q = q(w). The dispersion relation is physically measur- 
able and the same is true for the on-shell tensor I](w, q). 
For example, in the simplest case of transverse waves, the 
dispersion equation takes the form g^ = fcQl](a;,q) and 
the quantity E(w, q) can be referred to as the propaga- 
tion constant (index of refraction squared) of the Bloch 
mode^i^. 

We can seek approximate solutions to (|47|) by using 
the limit k — > 0. As was discussed in Sec. Ill B| the tensor 
I](a;,k) can be formally expanded in a non-gyrotropic 
medium according to (fT3)) . Wc substitute this expansion 
into (l47l) and obtain 



det [(k X (1 - /3)kx ) + fc^E(w, 0)] = 



(47) 



Equation ()47p is a valid approximation to the dispersion 
equation in the weak nonlocality regime and its solution 
yields the first nonvanishing solution to q (compared to 
the limit /i — > 0). Also, (|T7)) coincides with the disper- 
sion equation in a homogeneous medium with local pa- 
rameters e = S(a;,0) and /x = (1 — /3)~^. If we make 
this identification, we would arrive at the same homog- 
enization result as in the current-driven homogenization 
theory, except that we did not need to introduce the un- 
physical external current. However, this identification 
is not mathematical justified due to the reasons already 



discussed by us in Sec. IIIBI Here we reiterate these ar- 
guments in the somewhat new light of the Bloch-wave 
analysis. 

Firstly and most importantly, it can be easily seen that 
multiplication of (|47| by a scalar ^ does not alter the dis- 
persion equation or the value of q but doing so does alter 
the impedance J2°. Therefore, the above procedure is not 
expected to yield a meaningful correction to iF. This 
was illustrated above in Fig. 1111 In fact, S-parameter re- 
trieval predicts a much more accurate j2° while keeping 
approximately the same dispersion relation. This was 
illustrated in Fig. [71 And in general, it could not be 
reasonably expected that a theory that considers infinite 
media and disregards the physical boundary would pre- 
dict the impedance correctly. 

Second, the procedure described above is clearly inap- 
plicable outside of the weak nonlocality regime and, in 
particular, when ||/3|| ~ 1. In this region of parameters, 
introduction of the local permittivity and permeability 
tensors does not result in a correct dispersion relation, 
even approximately. But this is exactly the region of pa- 
rameters where current-driven homogenization predicts 
magnetic resonances. Consequently, this prediction is 
mathematically unjustified. 
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FIG. 15: Same as in Fig. [5] but for Example B. 



VII. DISCUSSION 

A. Current-driven excitation model and the theory 
of nonlocality 

The current-driven homogenization theory is deeply 
rooted in the theory of natural electromagnetic nonlo- 
cality (spatial dispersion) ^^1'^^ . The latter is, of course, a 
very successful theory, which has predicted and described 
theoretically such diverse phenomena as optical activity, 
additional waves, and anisotropy of crystals with cubic 
symmetry, etc. However, current-driven homogenization 
and, more generally, current-driven excitation model take 
certain analogies too far or apply them unscrupulously. 

The basic idea behind the current-driven excitation 
model can be traced to Ref. |3^. We translate the rele- 
vant text from the Russian edition of this book (Moscow, 
Nauka, 1965, p. 34), using only a slight change of nota- 
tions: 

"Generally, the arguments of the ten- 
sor I](a;,k) are mathematically-independent. 
This fact follows already from the definition 
(1.6) [equivalent to Eas. l48l49l below (authors' 
comment)] but can be at times not entirely 
obvious. This is so because, in optics, one 
encounters very frequently wave propagation 



FIG. 16: Example B. Real (top) and imaginary (bottom) 
parts of T as functions of kx/ko for h/\o = 0.66. EX - exact 
result, CD - current-driven homogenization. Only the range 
of kx/ko is shown for which T computed by either method is 
not negligibly small. 



in the absence of sources in the medium it- 
self, in which case k is a function of uj; for ex- 
ample, for normal homogeneous plane waves, 
k = {Ld/c)h{Ld, s)s. But if k = k(a;), then the 
spatial dispersion appears to be indistinguish- 
able from the frequency dispersion. This ob- 
servation raises a question [about the physical 
nature of spatial nonlocality (authors' com- 
ment)] and the answer to this question is the 
following. The tensor S](a;,k) is introduced 
for fields of the most general form, obtained 
when the sources Joxt(r) and /3ext(r) spatially 
overlap with the medium. Under these con- 
ditions, it is possible to create a field E with 
arbitrary and mathematically independent w 
and k (the Fourier components E(cj, k) is 
ultimately expressed in terms of Jext(w,k) 
and pext(w,k); see $2.1). From this, it fol- 
lows immediately that all problems involving 
wave propagation can be solved if S(w, k) is 
known." 

The last sentence in the above is only partially correct. 
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FIG. 17: Example B. Same as in Fig. [TB] but for R. 



In the case of natural nonlocality, S(w,k) is the spa- 
tial Fourier transform of the influence function (7(0;; r, r'), 
which appears in the nonlocal relation between D(a;,r) 
and E(a;,r), viz,: 



D(u;,r)== / cr(w; r,r')E(w,r')dV 
Jv 



(48) 



If both points r and r' are sufficiently far from the bound- 
ary of the medium, we can write ^(w; r, r') = /(w, r — r'). 
Then E(a;,k) is defined as the spatial Fourier transform 

of/(c^,r):' 



S(w,k)= f f{u},r)e'^''d^r 



(49) 



But this is insufficient to solve the boundary-value prob- 
lem for any finite shape. To that end, we would need 
to know how (t{lu; r, r') behaves when at least one of the 
points r and r' is close to the boundary. One can consider 
an approximation of the type 

a(c^;r,r') = 5(r)/(c., r - r')^(r') 

+ [1 - Sir)]6ir - r')[l - S'(r')] , (50) 

where S{r) is the shape function: it is equal to unity in- 
side the medium and to zero outside (in vacuum) . If (|50p 
holds, then the statement under consideration is correct: 
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FIG. 18: Example C. Same as in Fig. [S] but for Example C. 
Note that, in Example C, h is fixed while Ao varies. 



Maxwell's equations can be written in a closed form us- 
ing only the Fourier transform of /(w, r) and the Fourier 
transform of the shape function. Therefore, if E(cj, k) is 
known, then Maxwell's equations can be solved in a finite 
sample, at least in principle. We note that the familiar 
relation D(w,k) = I](w,k)E(aj,k) does not hold in this 
case and Maxwell's equations, written in the k-domain, 
contain an integral transform and cannot be solved by 
algebraic manipulation. This difficulty is known and ex- 
plained in $ 10 of Rcf. [38J but appears to be scarcely 
appreciated in the modern literature. But regardless of 
this difficulty, there is no reason to believe that ([SO]) is 
generally true. This approximation can be applicable, 
perhaps, if the nonlocal interaction between two points 
is transmitted only along the line of sight and if the body 
is convex. However, the first of these assumptions is dif- 
ficult to justify. 

The same analysis applies to current-driven homoge- 
nization of periodic composites. The knowledge of the 
function I](a;,k), as defined by ([5|) or by P5)) . allows one 
to find the law of dispersion but is insufficient to solve 
any boundary value problem. Therefore, this function is 
not an intrinsic physical characteristic of a composite. It 
is, rather, an auxiliary mathematical function, which ap- 
pears when a certain ansatz is substituted into Maxwell's 
equations written for an infinite periodic medium. 
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Now the difference between tlie classical theory of non- 
locality and the current-driven homogenization theory 
becomes apparent. In the former case, the real-space 
influence function (7(0;; r,r') is derived from first prin- 
ciples (e.g., from a microscopic theory) or introduced 
phenomenologically and then it completely characterizes 
the electromagnetic properties of a macroscopic object 
of any shape in the sense that it renders Maxwell's equa- 
tions closed. As discussed above, under some limited 
conditions, the dependence of (j{uj; r, r') on the two vari- 
ables r and r' can be simplified, e.g., as in ([50)1 . and 
then there exists a one-to-one correspondence between 
ct(w; r, r') and S(w, k). But in the case of current-driven 
homogenization theory, the low-pass filtering ([5]) does not 
follow from any first principle. Moreover, if ([5]) is ac- 
cepted as the fundamental definition of I](aj,k), there is 
no way to establish a one-to-one correspondence between 
the latter and the real-space function (j{uj; r, r'). As a re- 
sult, the knowledge of S(aj, k), thus defined, is insufficient 
to solve a boundary-value problem in any finite sample. 

Another obvious distinction between the two theories 
is that, for natural nonlocality, the infiuencc range (the 
characteristic value of |r — r'| for which cr(cj;r, r') is 
not negligibly small) is of the order of the atomic scale. 
Therefore, in the optical range, the condition of weak 
nonlocality is satisfied with extremely good precision and 
all the effects of nonlocality are, essentially, small pertur- 
bations. In the case of current-driven homogenization, 
the infiuencc range is h and the effects claimed as a re- 
sult of current-driven homogenization (e.g., /i w — 1) are 
dramatic and nonperturbative. 

So far, we have discussed the physical and mathemat- 
ical meaning of the function I](aj, k) as it is used both in 
the theory of natural nonlocality and in current-driven 
homogenization. The next important point to consider 
is the unjustified assumption of the current-driven ho- 
mogenization theory (and, more generally, of the current- 
driven excitation model) that the external current of the 
form ^ can, under some unspecified conditions, be cre- 
ated in the medium. This assumption also grows concep- 
tually from the above quote. In reality, "wave propaga- 
tion in the absence of sources in the medium itself" is en- 
countered in optics (and, more generally, in macroscopic 
electrodynamics) not just "very frequently," but always, 
without any known exceptions. Of course a medium can 
be optically active and emit some kind of radiation from 
its volume. However, in all such cases, the current inside 
the medium is subject to (linear or nonlinear) constitu- 
tive relations and cannot be created or controlled by an 
experimentalist at will. 

The above misinterpretation per se could have been 
of a purely academic or esoteric nature, but it has led 
to the unfortunate conclusion that a Fourier component 
E(cij, k), taken at an arbitrary value of k, is a physically 
realizable field and one can use it to compute certain ob- 
servable quantities that are quadratic in the field, such as 
the point-wise heating rate^^. This claim was addressed 
in our previous publicatioi*^. Here we note that E(a;, k) 



can be a term in an integral expansion of the real-space 
field E(a;,r). However, Plancherel's theorem notwith- 
standing, the expansion term E(a;,k) exp(ik • r) cannot 
by itself be used to compute any observable defined at a 
point r. 

Finally, another relevant misconception, which has 
been widely popularized in recent years^"— , is the 
proposition of equivalence of weak nonlocality of the di- 
electric response and nontrivial magnetic permeability. 
We have discussed this point in this article in much detail 
and have demonstrated that the equivalence exists only 
for the dispersion relation but not for the impedance of 
the medium. It can be argued that this is exactly what 
was meant by Landau and Lifshita^ since none of the rel- 
evant chapters consider boundary conditions in any form. 
The current-driven homogenization theory has taken this 
statement of equivalence out of its proper context and 
applied it to the problem of homogenization wherein the 
boundary conditions play the central role. 

It can be concluded that the classical theory of spatial 
dispersion is concerned primarily with certain physical 
effects such as rotation of the plane of polarization or 
appearance of additional waves, which are not present in 
the purely local regime but can be described as pertur- 
bations if small nonlocal corrections to the permittivity 
tensor are taken into account. The corrections are either 
introduced phenomenologically or computed using a mi- 
croscopic theory. The theory of spatial dispersion was 
never meant to be used for rigorous solution of bound- 
ary value problems and, therefore, the discussion of the 
dispersion equation sufficed in the vast majority of cases. 
For this reason, certain remarks appearing in the classi- 
cal texts on the subject, such as the now famous remark 
of Landau and Lifshitz regarding the equivalence of non- 
locality and magnetism, apply only to the dispersion re- 
lation. In the case of the current-driven homogenization 
theory, all these limitations have been disregarded. 



B. Homogenization by spatial Fourier transform 

Throughout the paper, we have emphasized the critical 
importance of taking the boundary effects into account in 
electromagnetic homogenization, particularly in the case 
of metamaterials whose lattice cell size typically consti- 
tutes an appreciable fraction of the vacuum wavelength. 
Consequently, theories relying entirely on the bulk be- 
havior of waves cannot be accurate; they may be capable 
of finding the effective index but not the impedance. 

Generally speaking, Fourier-based homogenization 
theories should be applied with extreme care because 
Fourier analysis makes it difficult to account for material 
interfaces, which break the discrete translational invari- 
ance of a periodic composite. 

It is feasible to devise a theory in which Fourier analy- 
sis is applied in the medium and in the empty space sep- 
arately. In this case, however, only natural Bloch modes 
will exist in the material, and no other values of k will 
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appear. The intuitive perception that a small localized 
source (e.g., a nano-antenna) embedded in a composite 
(e.g., in one of the empty voids) would generate within 
the material the whole spectrum of waves with all possi- 
ble real-valued k's is correct only in a very narrow tech- 
nical sense. In fact, within any area away from the small 
source, the waves with various k will interfere to produce 
the natural Bloch modes of the periodic structure. Only 
the latter are physically measurable. 



effective permittivity e rather than to an accurate predic- 
tion of /i. We note in passing that the correction to the 
magnetic permeability produced by the current-driven 
model is asymptotically 0{kh)'^, which is different from 
the 0{kh)^ asymptote that follows from S-parameter re- 
trieval. 



VIII. SUMMARY 



C. Current-driven homogenization 

As an illustration of the general principles stated 
above, we have critically analyzed the current-driven 
homogenization theory of Refs. |24||251 which does not 
account for the boundaries of the medium but derives 
the effective parameters from the behavior of waves in 
the bulk. In addition, this model relies on the use of 
physically-unrealizable sources inside the medium, with 
no justification as to why the results thus obtained should 
be experimentally relevant. 

The numerical results of Sec. |V] are therefore not sur- 
prising. They demonstrate that current driven homoge- 
nization does not yield accurate results in the range of 
parameters where it predicts nontrivial magnetic effects. 
In particular, the errors of the transmission and reflection 
coefficients T and R are of the same order of magnitude 
or, in some cases, much larger than the deviation of the 
magnetic permeability from unity, |j/i — 1||, where fi is 
the local permeability tensor predicted by current-driven 
homogenization. 

In most cases considered, current-driven homogeniza- 
tion does not provide a noticeable improvement in ac- 
curacy compared to the standard homogenization result 
([T6l). In the cases when such improvement can be ob- 
served, e.g., in Fig. [SI this is due to a correction in the 



This paper has four main conclusions. 

First, careful consideration of boundary conditions is 
required in all effective medium theories (EMTs). 

Second, all EMTs have an applicability range, and 
wherever a homogenization result is obtained, it is im- 
portant to verify that the parameters of the composite 
are within this range or, otherwise, validate the result 
with direct simulation. 

Third, there are many EMTs that yield the standard 
homogenization result in the limit ft, — >■ but different 
/i-dependent corrections to the former. Validating that 
these corrections are physically meaningful requires con- 
sideration of finite samples and cannot be done by inves- 
tigating an infinite periodic composite (this conclusion is 
closely related to the first one). 

Finally, current-driven homogenization theory fails to 
determine the impedance of the medium correctly and 
has, therefore, no useful predictive power, except in the 
limit /i — > 0, where it reproduces the standard homoge- 
nization result. 
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Appendix A: Determination of the coefficients m 
!J§a, s^b and Sh from the boundary conditions 



Upon substitution of the expressions ([20]) into (|2T|). vife 
obtain the foUovifing set of linear equations: 



[e^^Aa + e-^'f'-Ba 



A;, 



^e-"'" e'^"^a 



Kh 



1z 



Bb 



Qz 



Aa+B. 
An - B. 



Kb 



"^"Ah 



-"''"Bb 



'B, 



= 1, 
= 1, 

= 1, 
= 1 . 



This set is somewhat more complicated than what is en- 
countered in the ordinary theory of one-dimensional pho- 
tonic crystals. In the latter case, the matrix is the same 
but the right-hand side is zero. Correspondingly, the task 
is to find the value of kz (for a given kx) such that the 
equations have a nontrivial solution. This occurs when 
kz ~ Qz, where g^ is defined in Eq. (|24p . In this manner, 
the natural Bloch wave vector q of the medium is de- 
termined. For the case at hand, both k^ and k^ are free 
parameters but the right-hand side is nonzero. Therefore, 
the current-driven homogenization theory, essentially, re- 
places the problem of funding the natural Bloch mode of 
the medium by a mathematically unrelated problem of 
inverting the matrix in the above set of equations. 
The solution to the set stated above is given by 
where ^ is defined in (l23l) and 
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It can be seen that ^a is obtained from i/a (or vice versa) 
by changing the signs of the propagation constants k^ 
and Kb and of the related phases (pa and (pb (but not of 
6a and 6b), and similarly for the pair of coefficients ^4, 
3§b- Also, the coefficients are invariant under the the 
permutation of indices a ■H- b. 



Appendix B: Closed-form expressions for the EMPs 
obtained by current-driven homogenization 

In this Appendix we give the closed-form solution for 
the general case a ^ b. If a = 6, these expressions are 
significantly simplified. However, in the numerical codes 
used to produce the figures for this paper, we have used 
the general expressions given below. 

In addition to (fT5)) and (flB)) . we need to introduce the 
following notations: 

(i) The phase shifts (pa and <pb computed at k^ — are 
denoted by %pa and ipb'- 

i^a = Pa\k^=Q = Ka , Ipb = 4>b\k^=f) = hb ■ 

(ii) The refractive indices of each layer are denoted by 
ria = \feZ. and rib ~ ^Jlb- The branches of all square 
roots are defined by the condition tt < arg(-y/i) < 0. 

(iii) The dimensional size parameters for the lattice: 

Xa = k^a , Xb = kgh , X — Xa + Xb — koh . 

Note that ipa = riaXa, 'ipb = nbXb- 

(iv) The following symmetric combinations of the 
trigonometric functions: 

ai = 2(naSin-(/;a +?T-6sinV'b) , 
02 = 2 {uaCb sin -0a + UbCa siu Ipb) , 
as = Xa sin tpa + Xb sin ipb , 
a4 = xlua sin -0a + xlnb sin -ipb , 



and 

4i = Ua COS — sm — + Ub cos — sm — , 

^ ^b . i>a , i>a . Ipb 

6 = na COS — sm — + Ub cos — sm — , 

^3 = Ha sin Ipa COS 0b + Tif, sin 06 cos ^pa , 
^4 = nanfc-2-- — (?/;a sin -ipa COS 06 + i'b sin 0(, cos 0a) 
+ eaCfc (0b sin 0a cos ^pb + i^b sin 0h cos 0a) , 

<; / ^2 ■ 0a . 0b 

45 = (ea - Eb) sm -^ sm y , 



and 

r/ = ea(2 + cos0;,)sin^ -^ + 6^(2 + cos0a) sin^ — 



■nanb sin 0a sin 0b 



,,:. 2 0a „:„ 



C = Ea (40beb - 1pb<ia + SlpaUaUb) siu — siu 0f, 

+ efc (40aea - ijjaeb + SipbUarib) sin^ — sin 0a 
+ 8(6a-6b)^sin2^sin2|i 

0^?sin2%+0^^in^% 



and the following combinations of the functions in- 
troduced above: 

a = CatbiuA - Au'i) - a2 , 
r = 4(6^5 +66) . 

Note that by symmetry we mean here invariance 
with respect to permutation of indexes a O 6. 

Then the coscd-form expressions for the current-driven 
EMPs for the one-dimensional layered medium consid- 
ered in this paper can be stated as follows [we adduce 
the expressions for pi^x and pizz] fJ-xx and fXzz are given 
by dMl)]: 

<^yy = ^(eaeb)^^o , 

Pxx = yi-(efc - fa)^ [C - x'^PaPbiaebv] , 
6 



Pzz — --^{<^b - Ca) 



[naUb [2p + a;6(o' - aie||)] - Sai^s + t] , 
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and i^o is the determinant (|23p evaluated at k = 0: 
^0 = 1- cos(g^/i)|^^^o 

= 1 - cos 0a cos 06 - - 1 -^ H ) sin 0a sin 06 , 

2 V"b 'la/ 

and, finally, 

P = eaebS,2x'^{2+paPb2>o) 

"pa '4'b t 1 \ 

- naUbXsm — sin — (eae6X -I- 2e||^oj ■ 



Appendix C: S-parameter retrieval techniques used 
in this paper 

S-paramctcr retrieval is a well-established technique. 
However, the vast majority of papers that consider S- 
parameter retrieval arc focused on normal incidence only. 
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This limitation docs not allow one to access all tensor el- 
ements of the effective parameters. In this Appendix, we 
will remove this limitation and include off-normal inci- 
dences into consideration. 

Consider an anisotropic homogeneous slab charac- 
terized by the tensors e — diag(ej^, e^, e||) and /i ~ 
diag(^_L, /i^, /i||). The transmission and reflection coeffi- 
cients are given in Sec. IIVI It is convenient to rewrite the 
expressions given in that section in the following form: 



Here 
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p = cxp (iq^L) 



We refer the reader to Sec. IIVI for relevant nota- 
tions. Note the symmetries R{p, C) — —R{p, l/C") and 
T{p,C)=T{p,l/C). 

If T and R are known at some incidence angle (param- 
eterized by kx)-, so are C and p. The expressions for C 
and p in terms of T and R [inversion of (jCip ] is unique 
up to the branch of a square root and widely known: 
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where 
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In the case of low transmission (small |r|), one can use 
the approximate formulas 



C = 



C = 



1 + R 
^~R ' 
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p = 



(first branch) , 



1 -i?2 

1 — R^ 
p = (second branch) 



Here we have encountered the first instance of a branch 
ambiguity, but it can be easily resolved by using the con- 
dition of medium passivity, \p\ < 1. 

Other than the branch ambiguity mentioned, the above 
inversion formulas are well-posed and stable and allow to 
establish a one-to-one correspondence between the com- 
plex pairs {T,R) and {C,p). Given this result, we can 
assume that the functions C{kx) and p{kx) arc known as 
functions of k^ and will seek the effective tensors e and 
/i that fit these functions best at normal and closc-to- 
normal incidence. 

As is well known, the inverse problem posed above 
is ill-posed. We shall now describe three different ap- 
proaches to solving this problem. 



Let t : 
equation 



where 



Method 1 

kx/ko and x — k^L. Given the dispersion 
we can write 



C{t) 



Q{t) 



Pit) 
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Here n^ = e||/X|| is the squared index of refraction and 
recall that rj refers to /x for s-polarization and to e for 
p-polarization. 

Assuming C{t) and p{t) are known functions, we can 
express Q{t) by inverting each equation in (|C2[) . This 
will yield two solutions for Q{t): 



L(t) = ^^^7V^, Q2{t) = -\n[p{t)] 
C (t) IX 



At 



where A = 27r/.T and m an arbitrary integer. 

If wc could find such tensors e and /i that Qi{t) = Q2{t) 
for all values of t, then these parameters would describe 
transmission and reflection by the slab for all angles of 
incidence with perfect precision. This is possible only in 
the /i — > limit. Here we will require that the functions 
coincide at normal incidence and have the same second 
derivative with respect to t (the first derivative is identi- 
cally zero). To this end, it is convenient to introduce the 
functions 

Since these functions are expressed in terms of C(i), p{t) 
and known analytical functions (they do not contain any 
unknowns) , we can view them as directly measurable (or 
computable in terms of T and R) . Then the main equa- 
tion we wish to fit takes the form 

rji\F{t) = G(t) + Am . 

Here r^n and m (the branch index) are the unknowns. Of 
course, this equation has no solutions in general. But 
we can require that it holds to second order in t in the 
vicinity of i = 0. We can expand F{t) as 



F{t) = Fo + F2t^ , Fo = F{0) , F2 = lim 



F{t) - Fo 



and similarly for G{t). This results in a pair of algebraic 
equations 



V\\Fn 
V\\F2 



Go 

G2 



Am 



Even though this is a system of two linear equations with 
respect to two unknowns, it still cannot be solved because 
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m is, by definition, integer. We can, however, solve the 
first equation exactly (this will guarantee the correct T 
and R at normal incidence) and then minimize the norm 
of the second equation. This results in the following in- 
verse solution: 



m = Nint 



F0G2 — G0F2 



AF2 



At 



Fn 



(C3a) 
(C3b) 



where Nint[z] denotes nearest integer to the real part of 
z. Note that the second expression in the above (for 
7711) must use the value of m computed using the first 
equation. 

Once the quantities 7711 and m are known, we can also 
compute the squared refractive index according to 

n^ ^ (Go + Am)2 

Then, if ryu = /X|| (s-polarization), we can compute e|| = 
n-/jin. Otherwise, if r/ii = en. wc can compute /iii = 
nV?7||. 

This Method 1 docs not give access to e± and ij.±. To 
obtain these tensor elements, Method 3 must be used. 



incidence. We will fit the the index of refraction using 
second and forth normalized derivatives of p(t) and then 
use the function C{t) to find the impedance. 

We start by noting the following relations: 
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We can exclude the ratio ??||/77j^ from the above set of 
linear equations and solve for the index of refraction, viz. 



x{Fi/3F^ -1) ■ 
This gives the index of refraction in terms of the "mea- 
surables" F2 and F4. We can relate F2 and F^ to "mea- 
surements" of p(t) at some small but nonzero values of i, 
say, Ti and T2, as follows: 



„ t|6i - rf 62 

t2 - T 5 1 -^4 
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hi 



rl 



Method 2 



where 



This method combines Method 1 with the main idea 
of Ref. [S^. Namely, starting from some sufficiently low 
frequency a; or some sufhciently small /i, we gradually 
increase the relevant parameter (a; or K) and run apply 
Method 1 within this loop. However, after a few initial 
iterations of the loop, we change the rule according to 
which the branch index m is computed, namely, instead 
of using (|C3a[) . we chose the index m in such a way as to 
minimize the jump in the refractive index n. By "jump" 
we mean here In,- — n,_i I, where i is the iteration index. 



bh = -^ 



p{n) 
p{0) 



- 1 



fc = 1,2 



At this point we have found the index of refraction, n. 
We then compute 7711 and ri±^ from 



77|| = C(0)n , r]j_ 



.xCjO) 
F2 



Method 3 

This method is free from the branch ambiguity but 
does not guarantee exact fitting of T and R at normal 



In the second equation above, we have used rjn/rji^ = 
inF2/x^ as follows from (|C4ap . 

By considering both s- and p-polarizations, we can find 
all elements of the effective tensor. 



